stars <- function(x, digits=5, wide=F, align.stars=F, significance.level=.05) {
  if(significance.level == .05) {
    s1 <- .05
    s2 <- .01
    s3 <- .001
  } else {
    s1 <- .1
    s2 <- .05
    s3 <- .01
  }
  if(align.stars) i <- " & " else i <- ""
  cf = sprintf(paste0("%.",digits,"f%s%s"), round(x[1],digits), i,
             ifelse(x[4] < s3,"$^{***}$",
                    ifelse(x[4]< s2,"$^{**}$",
                           ifelse(x[4]< s1,"$^{*}$","")
                    )
             )
  )
  if(wide) cf <- c(sprintf(paste0("%.",digits,"f"), round(x[1],digits)),
                   ifelse(x[4] < s3,"$^{***}$",
                          ifelse(x[4]< s2,"$^{**}$",
                                 ifelse(x[4]< s1,"$^{*}$","")
                          )
  ))
  if(align.stars) {
    se = sprintf(paste0("\\multicolumn{2}{c}{(%.",digits,"f)}"), round(x[2],digits)) 
  } else {
    se = sprintf(paste0("(%.",digits,"f)"), round(x[2],digits)) 
  }
  if(wide) se = sprintf(paste0("%.",digits,"f"), round(x[2],digits))
  c(cf,se)
}

make.table <- function(cf,labels,f="",add.line=T, digits=5, 
                       wide=F,group.coefs=T, align.stars=F,
                       add.stats,add.stats.names, significance.level=.05,
                       wide.n=1) { # cf: list of fitted objects that have a (coef, se, t, p) $coefficient matrix 
  if(ncol(labels)==3) break.labels <- T else break.labels<-F
  nm <- length(cf)
  if(wide) {
    str <- paste0(paste(rep("%s",4*wide.n),collapse=" & ")," \\\\")
    str2 <- paste0(paste(rep("%s",2*wide.n),collapse=" & ")," \\\\")
  } else {
    str <- paste0(paste(rep("%s",nm+1),collapse=" & ")," \\\\")
    if(group.coefs) str2 <- paste0(str,"[.5em]") else str2 <- str
  }
  
  if(align.stars) {
    str0 <- paste0(paste(c("%s",rep("\\multicolumn{2}{c}{%s}",nm)),collapse=" & ")," \\\\")
  } else {
    str0 <- str
  }
  
  nv <- nrow(labels)
  for(i in 1:nv) {
    if(wide) est <- matrix(character(3*nm),nm,3) else est <- matrix(character(2*nm),nm,2)
    for(k in 1:nm) {
      if("coefficients" %in% names(cf[[k]])) this.cf <- cf[[k]]$coefficients else this.cf <- cf[[k]]
      if(ncol(this.cf) == 3) this.cf <- data.frame(Coef=this.cf[,1],SE=this.cf[,2], t=NA,p=this.cf[,3], row.names = rownames(this.cf))
      if(labels[i,1] %in% rownames(this.cf)) est[k,] <- stars(this.cf[labels[i,1],], digits=digits, wide=wide, align.stars = align.stars, significance.level=significance.level)
    } 
    if(!wide) {
      if(break.labels) {
        cat(do.call(sprintf, c(list(str), c(labels[i,2],est[,1]))),"\n", file=f, append=i>1) # 
        cat(do.call(sprintf, c(list(str2), c(paste0("\\;\\;",labels[i,3]),est[,2]))),"\n", file=f, append=T)
      } else {
        cat(do.call(sprintf, c(list(str), c(labels[i,2],est[,1]))),"\n", file=f, append=i>1) # 
        cat(do.call(sprintf, c(list(str2), c("",est[,2]))),"\n", file=f, append=T)
      }

    } else {
      if(nm>1) { cat("More than one model - wide option not possible."); return(NULL) }
      cat(do.call(sprintf, c(list(str), c(labels[i,2],est[,1], est[,2], est[,3]))),"\n", file=f, append=i>1) # 
    }
  }
  if(add.line) cat("\\hline \\\\[-1.8ex] \n", append=T, file=f)
  if(! (missing(add.stats) | missing(add.stats.names))) {
    for(k in 1:length(add.stats)) {
      if(!wide) cat(do.call(sprintf, c(list(str0), c(add.stats.names[k],add.stats[[k]]))),"\n", file=f, append=T)
      if(wide) cat(do.call(sprintf, c(list(str2), c(add.stats.names[k],add.stats[[k]]))),"\n", file=f, append=T)
    }
  }
} 

